function [beta2, T, T_J, F, part] = SURPartiesRobustLagvar_LogGMM(Y, X, mean_a, mean_b, RAND, n, k)


% First estimation
beta=inv(X'*X)*X'*Y;   % Get parameter estimates
U=Y-X*beta;  % Get residuals


% Second Estimation
P2=X*inv(X'*(kron(eye(n),ones(2,2)).*(U*U'))*X)*X'; %Second step projection matrix
beta2=inv(X'*P2*X)*X'*P2*Y; %Second step parameter estimates
u2=Y-X*beta2;  % Second step residuals
U2=(reshape(u2',2,n))'; % Reshape second residuals
clear U P2;

% Compute COVAR matrix
COVAR=inv((X'*X)*(inv(X'*(kron(eye(n),ones(2,2)).*(u2*u2'))*X))*(X'*X));

% Conduct Wald test on beta parameters
[F, T, T_J] = Wald_GMM(beta2, COVAR, k); 

% Estimate partial effects
bpartial=[beta2(2:k); beta2(k+2:2*k)]; % Create beta vector excluding constant
mean_la=[mean_a;mean_a;mean_a];    %create matrices for partial effects of binary variables
mean_ha=mean_la;
mean_ha(1,4)=1;
mean_ha(2,5)=1;
mean_ha(3,6)=1;
mean_lb=[mean_b;mean_b;mean_b];    
mean_hb=mean_lb;
mean_hb(1,4)=1;
mean_hb(2,5)=1;
mean_hb(3,6)=1;

[part] = partialLagvar_LogGMM(beta2, bpartial, mean_a, mean_b, mean_la, mean_lb, mean_ha, mean_hb, U2, RAND, n, k);